FINITE-ELEMENT DISCRETIZATION OF STATIC 
HAMILTON-JACOBI EQUATIONS BASED ON A LOCAL 
VARIATIONAL PRINCIPLE 
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Abstract. We propose a linear finite-element discretization of Dirichlct prob- 
lems for static Hamilton-Jacobi equations on unstructured triangulations. The 
discretization is based on simplified localized Dirichlet problems that are solved 
by a local variational principle. It generalizes several approaches known in the 
literature and allows for a simple and transparent convergence theory. In 
this paper the resulting system of nonlinear equations is solved by an adap- 
tive Gauss-Seidel iteration that is easily implemented and quite effective as a 
couple of numerical experiments show. 



1. Introduction 

With the advent [OS88] and success of level set methods and its many applica- 
tions }OF03i Set99] to areas ranging from computational physics to computer vision 
there has been considerable interest in numerical methods for solving Hamilton- 
Jacobi equations, dynamic and static. For problems with complex geometries, or for 
problems on manifolds, there is a demand for methods that work on unstructured 
meshes such as triangulations. 

Three main directions of constructing discretizations on unstructured meshes 
can be found in the literature. First, there are methods that lift ideas of finite- 
difference upwinding and Godunov schemes from hyperbolic conservation laws to 
Hamilton-Jacobi equations (whose solutions are, at least in ID, integrals of solutions 
to conservation laws), see, e.g., [6898] • Second, there are finite-element methods 
that are based on a weak formulation of the semi-linear second order equation 
obtained by adding a small amount of factual viscosity, see, e.g., [LYC03] . And 
third, there are methods that utilize the connection of Hamilton-Jacobi equations 
(via Bellman's principle) to optimal control problems, see, e.g., }SV03] . 

In this paper we propose a discretization that bears similarities with the last 
approach. We implicitly construct a linear finite-element solution by requiring that 
it solves locally a simplified equation with (local) boundary conditions given by the 
finite-element function itself. The simplified local equation is then solved by a local 
variational principle, the Hopf-Lax formula. 

This simple discretization is interesting in various respects. First, we will show 
that it generalizes quite a few approaches known in the literature. Second, it allows 
for an extremely simple, self-contained convergence theory. In fact, the only results 
of the general theory that we rely on are a uniqueness theorem for viscosity solutions 
and the theorem of Arzela-Ascoli. Existence of viscosity solutions will be shown in 
passing by the convergence of the finite-element discretization. 

By construction the discretization inherits structural properties of the viscosity 
solution of the Hamilton-Jacobi equation such as a comparison principle. For each 
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property we will carefully trace the specific assumptions on the Hamiltonian and the 
boundary data that are needed for proofs in the continuous and the discrete case. 
In particular, it is known [BCD97i |Lio82] that the existence of a viscosity solution 
of the Dirichlet problem necessitates a compatibility condition on the boundary 
data, which is basically a restrictive Lipschitz bound. However, this necessary 
condition gets barely any mention in the literature on numerical methods, even in 
the formulation of convergence results, e.g., }SV03i Thm. 7.7]. 

Moreover, we propose in this paper a likewise simple iterative method for solv- 
ing the resulting nonlinear system of equations, namely the adaptive Gauss-Seidel 
iteration, which is easily implemented and, at least experimentally, quite effective. 

The paper is organized as follows. In |2] we recall the concept, existence, and 
uniqueness of viscosity solutions of Dirichlet problems for certain Hamilton-Jacobi 
equations. In fj3] we introduce the support function of the zero-level set of the 
Hamiltonian which plays a major role in the definition of the discretization. In 
we define the linear finite-element solution based on a local variational principle. 
The existence, uniqueness, and uniform Lipschitz continuity of the finite-element 
solutions are subject of §[5l A suitable concept of consistency is introduced in 
§6] and the convergence of the discrete solutions is proved. In JT] we apply the 
finite-element discretization to a class of generalized eikonal equations in 2D and 
obtain, by a simple geometric argument, a closed formula for the local discrete 
equation. In §8] we discuss the adaptive Gauss-Seidel iteration that we propose 
for an easily implemented and quite effective solution of the nonlinear system of 
equations. Finally, in §9] we study two numerical experiments and compare the 
proposed method to the ordered upwind method (OUM) recently published by 
Sethian and Vladimirsky [S V03] . 



2. Existence and Uniqueness of Viscosity Solutions 

In this section we shortly review the existence and uniqueness theory for the 
Dirichlet problem of a Hamilton-Jacobi equation, 

H{x,Du{x)) = 0, xen, u\on=g, (1) 

where throughout the paper we will assume that C K'' is a bounded Lipschitz 
domain. For convex Hamiltonians H a sufficiently general set of assumptions is 
(see pb82l §5.3]): 

(HI) (Continuity) H eC{TixR'^). 

(H2) (Convexity) p '—^ H{x,p) is convex for all x £ VL. 

(H3) (Coercivity) H{x,p) — s- oo as ||p|| — > oo, uniformly in x Q. Equivalently, 
by assumptions (HI) and (H2), there are positive constants a, (3 with 

H{x,p) a\\p\\ ~ f3, xen,peR'^. 

(H4) ( Compatibility of the Hamiltonian) H(x, 0) ^ for all x Q. 
Existence of a solution to (J]) requires a further condition on the boundary data: 

(H5) ( Compatibility of Dirichlet data) g(x) — g{y) ^ S{x, y) for all x,y dfl. 

Here, S denotes the optical distance defined, under the assumptions (H1)-(H4), by 

6{x,y)^miy\{at),-e{t))dt : ^ G C°^'([0, 1],!!), ^(0) = C(l) = 2/} 

where p(x, q) = max (p, q) . (2) 

H(x,p)=0 
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In fact, S qualifies as a distance by the fairly obvious properties 

6{x,x)=0, ^ S{x, z) ^ d{x,y) + S{y, z), x,y,z£Q.. (3) 

If H is symmetric with respect to p, that is, H{x,p) = H{x, ~p), then S defines a 
pseudometric on fl. 

Let us recall the concept of viscosity solutions [CEL84] for the first order equation 

H{x,Du{x)) ^0, xen. (4) 

A function u € C°'^(f2) is a viscosity suhsolution (super solution) of (3| if all 
V e C^(yi) with u — V attaining a local maximum (minimum) at some Xq & Q. 
yield 

H{xQ,Dv(xo)) i^Q (^0). 
Now, a viscosity solution is simultaneously a viscosity sub- and supersolution. Note 
that by Rademacher's theorem on the differentiability of Lipschitz continuous func- 
tions a viscosity solution satisfies (EJ) pointwise almost everywhere. 

Theorem 1 (P.-L. Lions tLio82i Thm. 5.3]). Assume (H1)-(H4). The Dirichlet 
problem (J\) has a viscosity solution u if and only if the boundary condition satisfies 
the compatibility condition (H5). A specific viscosity solution is then given by the 
Hopf-Lax formula 

u{x) = inf (g{y) + S{x, y)) . (5) 

While this theorem will only serve as a motivation for the finite-element dis- 
cretization in §5l we will obtain the existence of viscosity solutions under somewhat 
more restrictive assumptions as a spin-off of the convergence result, Theorem \Tll 

Uniqueness requires a compatibility condition on the Hamiltonian that is slightly 
stronger than (114): 

(H4') H{x,0) < for all x € n. 

In fact, uniqueness of the viscosity solution given by (E) is then a simple corollary 
of the following comparison principle. 

Theorem 2 (H. Ishii [Ish87] ). Assume (H1)-(H3) and (H4'). Letu,v be viscosity 
sub- and supersolutions of 0), respectively. If u ^ v on dCl then u on CI. 



3. The Support Function of the Zero-Level Set 

In the definition (2) of the optical distance S the support function (see }Roc70i 
p. 28]) 

p{x,q)= max {p,q) = sup {p,q), xen,qeM.'^, (6) 

H{x,p)=0 H{x,p)!ia 

of the zero-level set of H made an appearance. It is a well-defined real-valued 
function, since by (H3) the zero-level set is compact and by (114) non-empty. The 
second equality in (6J follows from the convexity (112). 

Since the discretization that we propose in the next section will be based on this 
support function p we collect its most important properties. 

Lemma 3. Assume (H1)-(H4). Then p : fix S.^ — > M is upper semicontinuous in 
the first argument, positively homogeneous convex in the second, that is, 

p{x,qi + q2) ^ p{x,qi) + p{x,q2), p{x,tq) ^ tp{x,q), t ^ 0, 

for X E n and q, qi,q2 G K'^. Let p* — (3/a with a, (3 from (H^). Then 

0^p{x,q)i^p*\\q\\, xeH, geM'^. 
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Assume additionally (H4-')- Then pIqxW continuous and 

p{x,q) > 0, xefl, O^qeW^. 

Proof. Being defined as the pointwise supremuni of linear functions the function 
q p{x, q) is a convex, positively homogeneous and, by (H4), nonnegative function 
for fixed x £ il. Assumption (H3) yields for H{x,p) ^ the bound ||p|| ^ f3/a = p* , 
which readily implies the upper bound on p. 

To prove the upper semicontinuity let x„ — + xq G O and q G M.^. We extract a 
subsequence Xn' such that 

p{xn',q) = limsupp(x„,g), 

n — *oo 

where p„' is a maximizing argument with H{xn' ,Pn') = 0. Because of the bound 
llPn'll ^ P* we can assume without loss of generality that p„' — s- pq. We obtain 
H(xq,po) = and therefore 

limsupp(a::„,g) = {po,q) p{xo,q). 

n — *oo 

From now on, we assume (H4'). Let a; G fi. Since H{x,0) < there is (5 > 
such that H{x,p) ^ for ||p|| < S. Thus, for g 9^ 

p{x, q) ^ max (p, q) ~ 5 \\q\\ > 0. 

\\p\\ 

Finally, to prove the lower semicontinuity let a;„ — > a;o G f2 and q G . There 
is a maximizing po G 1^'* with p{xo,q) — {po,q) and H{xo,pq) = 0. We extract 
a subsequence such that p[xn',q) ^ liminf„_^oo p(a;ri, <?) and, below, construct a 
sequence p„/ — > po with H(xn' ,Pn') ^ 0. With it in hand we conclude 

liminf p(a;„,g) = lim p{xn-,q) ^ lim (p„',g) = (po,g) = p{xo,q)- 

n — ^00 n' — ^00 n' — ^cxd 

There is no loss of generality in assuming that either always H{xn',Po) ^ or 
always H{xn',Po) > 0. In the first case we simply take p„/ = po- In the second 
case, since H{xn' , 0) < 0, there is a A„' G (0, 1) with H{xn', A„'Po) = and we put 
Pn' = An'PQ. We can assume that A„' — > Aq G [0, 1]. Taking limits in 

= H{Xn>,Xn'Po) S% (1 - Xn')H {Xn',0) + \n'H{Xn' ,Po) 

yields 0^(1 — \q)H{xq, 0) which, by (H4'), implies Aq = 1 and p„' — po. □ 

Note that if (H4') holds and H is symmetric with respect to p, then q ^ p{x, q) 
defines a norm on R'^ for all a; G ft. 

Lemma 4. Assume (H1)-(H4-) and that the segment joining the points y,z £ Q 
belongs to fl. Let p^, ^ be a constant such that p{x, q) ^ x G Vl and q G M''. 

Then, with the constant p* defined in Lemma [H 

p* \\y- z\\ ^ 5{y,z) p*\\y- z\\. 

Lf H{x,p) does not depend on x, then p{x,q) = p{q) does not depend on x either 
and 

5iy,z) = p{y - z). 
Proof. The optical distance S{y, z) is bounded by the expressions 

Po - inf 1^ U'{t)\\ dt : ^ e C"^i([0,l],n) such that C(0) = y, ^(1) = 

with pq = p^f for the lower bound and po — p* for the upper bound. The infimum 
is nothing but the minimal length of a path joining the point y and z within fl. By 
the assumption on y and z this minimum is realized by the segment joining them. 
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If H , and hence p, does not depend on x G fJ, wc obtain by Jensen's inequality 
for ^ e C'HfO, 1], f^) with ^(0) = y and ^(1) = z that 

p(-r(t)) rft ^ P dt^ - p(y - z). 

The lower bound is attained for the segment joining y and z yielding the assertion 
5{y, z) = p{y — z) (see also [Lio82l Remark 5.7]). □ 

Example. An important class^ of Hamiltonians satisfying (H1)-(H3) and (H4') is 
given by 

H{x,p)^F{x,p)-l 
where F G C{Ti x R^) is assumed to be positively homogeneous convex in p with 
the bounds 

< F{x,p) F*, xen,\\p\\^l. (7) 
Duality theory of nonnegative positively homogeneous convex function (gauges) 
[Roc70l §15] teaches that the support function p of the zero-level set of H is the 
polar of F, that is, 

p(x, q) — max — r, r [x,pj — max. ■ 



P#o F{x,p)' ' q^o p{x,q)' 

Hence, for the (particular) Hamilton-Jacobi-Bellman equation jSV03i Eq. (22)] 
with 

H{x,p) ^ max{p,-q)f{x,q) - I, (8) 

Il9ll = l 

where / is continuous with bounds < /* ^ f{x,q) ^ /*, x G \\q\\ — 1, we 
immediately read off that 

P{x,q)^ "'^''i II, , xen,qeR''. (9) 

f[x,-q/\\q\\) 



4. The Finite-Element Discretization 

Linear Finite Elements. Let us shortly recall the notion of linear finite-elements. 
For a sequence ft, ^ we consider a family E/j of shape-regular simplicial triangula- 
tions of (the now polytopal domain) f7 C R'^. We denote the diameter of a (closed) 
simplex a € T,h by fti(cr) and the minimal height of a vertex in a by hQ{a). We 
assume 

h — max hi(a) 
and measure the shape-regularity by a uniform bound 

ho{a) 

where we call 9 the regularity constant of the family of triangulations. 

The space of linear finite elements on E/i, that is, continuous functions that are 
affine if restricted to a simplex cr G S;i, is denoted by Vh and 

h ■■ cm Vh 

is the corresponding nodal interpolation operator. We endow Vh with the maximum 
norm, that is, convergence in Vh is the uniform convergence of the finite-element 
functions. 

The set of nodal points (vertices) of the triangulation Y^h that belong to fl, il, 
dQ are denoted by flh, f^/i, dflh, respectively. Note that a finite-element function 



Which essentially covers the general case as wo will see in Footnote [2] in |6J 
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Figure 1. The neighborhood iJh{xh) of Xh € ^h, that is, the collection 
of all simplices (like the one shaded in dark) that have Xh as a vertex. 



Uh G Vh is uniquely determined by its nodal values, that is, the values Uh{xh) for 
all e 0?i. 

For an interior nodal point xt G we consider the simplicial neighborhood 
'^h{xh)i that is, the interior of the miion of all simplices in Y,h that have Xh as a 
vertex (see Figure QJ- 

The Idea. The finite-element discretization which we propose is motivated by the 
idea of local solutions: 

At Xh € the finite-element solution Uh G Vh takes the value 
u*f^{xh) of the exaci viscosity solution mJ^ e C^'^{ujh{xh)) that solves 
a simplified Hamilton- Jacobi equation on LUh{xh) subject to the 
boundary conditions = Uh\auj^{xh)- 

A good candidate for such a simplification of the Hamilton-Jacobi equation {IJ is 
obtained by freezing locally the dependence of H on its first variable. This way 
uf^ G C'^'^{LL!h{xh)) is obtained as the viscosity solution of the local Dirichlet problem 

H{xh,Du*h{x)) =0 on uj{xh), u*h\9^,^(^^^) ^ Uh\duj^{x^)- (10) 

This simplification is particularly suitable, because for [y,z] C tOhixh) the optical 
distance Sx^{y, z) of the local equation (JM is, by Lemma 3^ 

Sxh{y,z) = p{xh,y- z), 

where p{xh,-) is the support function of the zero-level set of the convex function 
H{xh,-) as defined in (6}. Theorem Q] tells us that if the local Dirichlet problem 
( [To] ) is solvable, then the value u'^(xh) is given by the Hopf-Lax formula 

u*h{xh) = min (uh{y)+p{xh,Xh-y)), 

that is, by a simple local variational principle. We note that, under the assumptions 
(H1)-(H4), the Hopf-Lax formula is well-defined independently of the compatibility 
of the boundary data of the local Dirichlet problem ( [TO] ). Anyway, we base the 
finite-element discretization on this formula and the convergence will be proved 
later without the interpretation of as a local solution. 

The Discretization. We define a function A^ ■ Vh Vh, called the Hopf-Lax update 
function, by 

{min (uh{y) + p{xh,Xh-y)), Xh e ^h, 

Uh{xh), Xh G dVlh. 
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The finite-element solution Uh £ Vh that discretizes the Dirichlct problem (jl]) is 
now implicitly defined by the fixed-point equation 

Uh ^ AhUh, w/i|anh = fflaofc- (H) 

As in the continuous case, we call Uh £ Vh finite-element subsolution (supcrsolu- 
tion) if Uh s% AhUh [uh ^ AhUh). 

We remark that the evaluation of (AhUh){xh) at an interior nodal point Xh G 
can be calculated by a finite collection of (d — l)-dimensional convex optimization 
problems. This follows from the representation (see Figure Q]) 

{AhUh){xh) = min min \ Uh{y) + p{xh, Xh ~~ y) ■ 

y G the {d — l)-dimensional face of a opposite to Xh| (12) 

and the observation that Uh is affine, and hence Uh + pixh, Xh — •) convex, on a. In 
§[7] we will study an important class of examples for which these {d— l)-dimensional 
convex optimization problems allow for a particularly simple solution. In general, 
however, one would have to use suitable iterative numerical methods to solve them. 

Remark. For the particular Hamilton-Jacobi-Bellman equation |8]) the finite-ele- 
ment discretization pT| ) is equivalent to various grid-based methods that are ob- 
tained from linear interpolation of the grid values and a direct local application of 
Bellman's dynamic programming principle. See, e.g., [Tsi95l Eq. (2.3)] and [SV03i 
Eq. (25)] as well as the references given therein. 



5. Existence and Uniqueness of the Finite-Element Solution 

The existence of a finite-element solution as implicitly defined by (TTj ) is based 
on two simple properties of the Hopf-Lax update function Ah- 

Lemma 5. Assume (H1)~(H4). Let Uh,Vh G Vh- 

(1) Ah is monotone, that is, Uh ^ Vh implies AhUh ^ AhVh- 

(2) Ah is nonexpanding, that is, \\AhUh - AhVh\\^ ^ \\uh - Vh\\^. 

Proof. The first property is an immediate consequence of the definition of Ah. To 
prove the second, let the maximum be attained at a nodal point Xh G ^k, without 
loss of generality \\AhUh - A;,w,,||^ {AhUh){xh) - [AhVh){xh)- If Xh e dVlh there 
is nothing to show; so we can assume Xh £ ^h- Let y^^ G dcuhixh) be such that 

{AhVh)ixh) = Vh{y*) + p{xh,xh - y*). (13) 

Hence 

{AhUh){xh) - {AhVh)ixh) 

^ (uhiy*) + p{xh, Xh - y*)j - (vhiy*) + p{xh,xh ~y*)j ^ \\uh -~Vh\\^, 
which proves the assertion. □ 

Theorem 6. Assume (H1)~(H4) and g : 9f2 —^ M. Then the finite- element dis- 
cretization pi] ) has a solution Uh G Vh. If u1 G Vh is such that u1\qq^ = glon^ 
and AhU^h ^ '"h' then the fixed point iteration 

ul+^=Ahul, n = 0,1,2,..., 

converges monotonously to a solution of pTl) . 
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Proof. An initial iterate e Vh with A^m^ ^ is given by 
ulldn^ = g\dn^, Uh\n,, = mm g{x). 

Inductively the monotonicity of A/j implies uJJ^^ = A/juJJ ^ uJJ. Hence, the mono- 
tone convergence of the sequence follows if we establish a uniform bound on the 
iterates. Since such a bound is trivial for the boundary we consider for a given 
Xft S fl/i a shortest path Xh = . . . ,x™ of nodal points that connects xu along 
edges of the triangulation with the boundary: S d^h- There is a bound L on 
the length of such a path which depends on the triangulation but not on Xh- With 
p* as defined in Lemma [3] we get 

m — 1 

<{xh) ^ max g{y) + ^ p{xl, ^ - x\^^) max g{y) + p*L. (14) 

4=0 

Thus, — *■ u;,. G V/i for some Uh € V/,., which by continuity must be a fixed point 
of A;,. □ 

Like in the continuous case, uniqueness of the finite-element solution requires 
the sharper condition {H4') and is a simple corollary of the following discrete com- 
parison principle. Thus, the finite-element discretization is a monotone scheme. 

Theorem 7. Assume (H1)-(H3) and (H4')- Letuh,Vh G Vh be finite- element suh- 
and super solutions, respectively. If Uh ^ Vh on dflh then Uh ^ Vh on Vl. 

Proof. Let be A/j = — G Vh. Note that the maximum of A/j will be attained 
in a nodal point. We will show that the existence of Xh G Oh with /S.h{xh) = 
max^gQ A/i(x) = 5 > yields a contradiction. To this end we choose such a 
maximizing Xh with minimal value of Vh{xh). With j/* G duJh{xh) as in ( IT31 ) we get 

S = Uh{xh) - Vh{xh) < {KhUh){xh) - {KhVh){xh) ^ Uh{y*) - Vh{y*). 

On the boundary of the face that contains in its relative interior there is, by the 
maximality of 5, a point x\ G 0.h such that Aft(x^) — S and Vhix'^) ^ Vh{y*). By 
Lemma [3] we have p(xh,Xh — y*) > and obtain 

Vh{x*h) < Vh{y*) = {KhVh){xh) - p{xh,xh - y*) < {KhVh){xh) ^ Vh{xh) 
in contradiction to the minimality of Vh{xh). □ 

In the discrete case, up to now, we did not impose a compatibility condition 
on the boundary data such as (HS). This will change in the discussion of a third 
important property of the finite-element solutions needed for convergence, that is, 
uniform Lipschitz continuity. Based on the constant ^ of Lemma 0] and the 
regularity constant ^ 1 of the family of triangulation we consider the condition 

(H5') g{x) - g{y) <c: ^ - y|| for aU a;, y G 9^!. 

By Lemma 3] this condition is actually stronger than (115). Note that the homoge- 
neous Dirichlet condition g = always satisfies (HS'). 

Theorem 8. Assume (H1)-(H3), (H4' ) and (H5'). The unique finite-element 
solution Uh G Vh of (11\ ) satisfies the uniform Lipschitz condition 

\uhix) -Uh{y)\ i^cndd- p* ■ \\x - y\\ , a;,?/ G H, 

and the uniform bound 

||tih||oo < max \g{x)\ + cnOd ■ p* ■ diam(ri). 

Here, 6 denotes the regularity constant of the family of triangulations, p* is the 
constant defined in Lemma [B\ and cq > is a constant depending only onfl. If fl 
is convex, we can choose cn — 1. 
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Proof. The uniform bound on is a simple consequence of the Lipschitz con- 

dition. The proof of the Lipschitz condition proceeds in three steps, imposing less 
and less restrictions on the possible choices oi x,y G Q. 
Step 1. For neighboring nodal points Xh,yh S we prove 

\uh{xh) ~ Uh{yh)\ ^ P* ■ \\xh - Vh\\ ■ 

Since p* ^ p* ^ p*/0 this is, by (H5'), obviously true for Xh^Vh G dVth- If Xh G 
we have yt G dujh{xh) and hence 

Uh{xh) = {KiiUh){xii) ^ Uh(yh) + p{xh,xh - yh) ^ Uh{yh) + P* \\xh - yh\\ ■ 

If yti € we can change the roles of x^ and yh and the Lipschitz bound follows. 

Assume on the other hand that yh G dilh- There is a minimizing G dujh{xh) 
such that 

Uh{xh) = {^hUh){xh) = Uh{y*) + p{xh,xh - y*) > Uh{y*), 

where the last inequality follows from Lemma \3\ The boundary of the face that 
contains ?/* in its relative interior has a point x\ G flh with Uhix\) ^ Uh{y*) < 
Uh{xh)- By the definition of and we obtain 

p{xh,Xh - y*) ^ p4xh - y*\\ ^ yll^'i ~ ^/ill- 
Continuing this construction we obtain a sequence Xh = x^,x\,. .. , a;™ of nodal 
points with strictly decreasing u^-values that necessarily reaches the boundary at 
some index to: G dilh- Thus, by construction and (H5'), 

m — 1 

Uh(xh)^g{x]:) + ^Y.\\'^h-<^'\\ 

i=0 




^ uhiyh) - y ||2;/i - yh\\ > Uhiyh) - p* \\xh - yh\\ , 
which concludes the proof of Step 1. 

Step 2. Let cr G be a simplex of the triangulation. For x,y £ a we prove that 

\uh{x) ~ Uh{y)\ s^9d- p* ■ \\x - y\\ . 
By an afhne transformation x Bx+b we map the standard c?-dimensional simplex 

a = {xe M^o : £i + . . . + s; 1} 

onto a. The puUback of Uh\a under the transformation will be denoted u. By Step 1 
we can estimate the length of the (constant) gradient of Uh\a- by 

\\Duh\a\\ ^ \\B-^\\ \\Du\\ ^ \\B-^\\Vdp* ■ hi{a). 

Now, is the largest ratio of the length of a segment in a to the length of 

its image in a. Without loss of generality such a segment can be assumed to join 
a vertex with the opposite boundary face. Thus ^ \/d/hQ{a) and, by the 

shape-regularity assumption, that is, /ii((T)//io(cr) ^ 0, we get 

\\Duh\a\\ i^0d-p* 

and hence the assertion of Step 2. 
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Step 3. For x,y € VL there is a Lipschitz path 7 e C°'^([0, 1], fi) joining x and y 
such that (see [Alt99l p. 304]) 

h'lloo «^ cn\\x-y\\. 

For convex the path 7 can be chosen as the segment joining x and j/, which yields 
CQ = 1. Now, let = < ^1 < • • ■ < = 1 be a subdivision of [0, 1] such that 
7(ii_i) and ^{ti) are elements of a common simplex. By Step 2 we obtain 

■m— 1 m— 1 

- ^ ^ |u,,(7(i,)) - Uh{l{U+i))\ Od ■ p* Ihi^i) ^ 7(^.+i)ll 

m — 1 

(^cnOd- p* ■ ||a;-y|| ^ ^cnOd- p* ■ ||x-y||, 

1=0 

which concludes the proof of the asserted Lipschitz bound. □ 



6. Convergence of the Finite-Element Discretization 

The argument will be simplified if we consider a modified Hamiltonian H for 
which the corresponding Hamilton- Jacobi equation possesses the same viscosity 
solutions as the original one. 

Lemma 9. Assume (H1)-(H4). Let x G Cl and p E M''. For the modified Hamil- 
tonian 

H{x,p) = max {{p,q) - p{x,q)) 

Il9ll = l 

we get that H{x,p) ^ (H{x,p) ^ 0) implies H{x,p) ^ (H{x,p) ^ 0). 

Proof. First assume H{x,p) > 0. There is a hyperplane that separates p strongly 
from the compact and convex level set {p : H{x,p) ^ 0} (see }Roc70i Cor. 11.4.2]). 
That is, there is a vector q € M.'^, \\q\\ = 1, with p{x,q) < {p,q). Hence 

H{x, p) ^ (p, (?) - p{x, q) > 0. 

Now assume H{x,p) < 0. There is e > such that H{x,p + 5p) < for \\6p\\ ^ e. 
Hence 

{P, q) - P{x, q) < {P, q) - , max {p + 5p,q) = -e\\q\\ , 
Taking the supremum over all q with \\q\\ = 1 yields H{x,p) ^ — e < 0. □ 

In particular, each viscosity subsolution (supcrsolution) of the thus modified 
Hamilton-Jacobi equation H{x, Du{x)) = 0, a; S i7, is also a viscosity subsolution 
(supcrsolution) of the original one H{x, Du{x)) = 0, a; € 51.^ 

Loosely speaking, in the framework of viscosity solutions the notion of consis- 
tency of a discretization means that a smooth function is already a subsolution 
(supcrsolution) of the differential equation if it is a subsolution (supcrsolution) of 
the discrete scheme. The precise statement is given in the next theorem. 



^Under the additional assumption {H4') the same holds true, since then p{x, q) > for x £ fl, 
if we consider the modified Hamilton— Jacobi equation with the Hamiltonian 

H(x, p) = max 1. 

\\q\\^l p{x,q) 

Hence, we see that the example at the end of |2] in fact covers the general case. 
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Theorem 10. Assume (H1)-(H3) and (H4'). Letv G Cg°{Q), x eQ, and Xh G flh 
be a sequence of nodal points that converges to x as h 0. Then 

v{xh) ^ {hhlhv){xh) for all h =i> H{x, Dv(x)) 0, 

v{xh) ^ {Ahlhv){xh) for all h =i> H{x,Dv(x)) ^ 0, 

where Ih ■ C'{fl) — > Vh denotes the nodal interpolation operator. 

Proof. Since v is smooth we can approximate the directional derivatives of v in Xh 
by first order differences as foUows 

Dvixh),jp^^)+Oih), yediJ^ix,,). (15) 

\\xh-y\\/ 

Now, let v{xh) ^ {Ahlhv){xh) for all h of the sequence, that is, 

v{xh) - (Ihv){y) - p{xh, Xh - y) 0, y ^ dujh{xh)- 

After division by ||a;;i — y\\ we get, by (15), a constant c > such that 

{Dv{xh),q) - p{xh,q) ^ ch, \\q\\ = 1. 

If we pass to the limit h ^ (note the continuity of p at a; G as stated in 
Lemma [3) and take thereafter the maximum over all \\q\\ ~ 1, we obtain 

H{x, Dv{x)) — max {{Dv(x), q) — p{x, q)) ^ 0. 

Il9ll = l 

From Lemma [9] we infer the assertion H{x, Dv{x)) ^ 0. 

On the other hand, let v{xii) ^ [Khlhv){xh) for all h of the sequence, that is, 

v{xh) - {Ihv){yh) - p{xh,xh ~yh) ^0 

for some y^ G dujhixh)- After devision by \\xh — yh\\ we get, by (15) , a constant 
c > such that 

(Dv{xh),qh) " p{xh,qh) ^ -ch, qt = [xh - yh)/\\xh - yh\\- 

By compactness, we can assume that qh — > q* with \\q^\\ = 1. Passing to the limit 
/i — > we thus obtain 

H{x, Dv{x)) ^ {Dv{x), q,) - p{x, q,) ^ 0, 

from which we infer the assertion H{x, Dv{x)) ^ by Lemma [91 □ 

Now we have all the tools in hand to prove the convergence of the finite-element 
discretization. 

Theorem 11. Assume (H1)~(H3) , (H4'), and (H5' ). Then, as h ^ 0, the sequence 
of unique finite- element solutions Uh G Vh defined by 

Uh = AhUh, Uhldflh = gldQh^ 
converges uniformly to the unique viscosity solution u of the Dirichlet problem 

H{x, Du{x)) = 0, u|af2 = g. 

Proof. Theorems [6] and [7] show the existence and uniqueness of the finite-element 
solutions Uh G Vh- Theorem [8] shows that Uh G Vh is a uniform bounded sequence 
of uniform Lipschitz continuous functions. By the theorem of Arzela-Ascoli there 
is a subsequence {uh') that converges uniformly to a function u G C°^^(ri). Because 
of (H5') and = fflaOhj this limit satisfies the boundary condition u\qq, = g. 

To show that u is a viscosity suhsolution of H{x, Du{x)) = let ?; G C§°{^1) 
and xo G such that u — v attains a local maximum in xq. By adding a quadratic 
parabola to v if necessary, we may assume that it is in fact a strict local maximum 
(see [Eva98 i p. 542]). Extracting a further subsequence of h' if necessary, there is, 
by uniform convergence and the monotonicity of the nodal interpolation operator 



Vhjxh) - {hv){y) 

\\xh - y\\ 
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Ih ■ C(ri) ^ Vh, a. sequence of nodal points Xh' S flh' such that Xh' Xq and (see 
the argument given in [Eva98, p. 541]) 

{uh> - v){xh') > (uh, - Ih'v){y), y E dujh'{xh')- 

Now let 2/* e du)h'{xh') be a minimizing argument such that 

{Kh'Ih'v){xh') = {Ih'v){y^) + p{xh',Xh' -y*)- 

Then it holds that 

[uf - v){xh') > Uh'(y*) + p{xh',xh' - y*) - {Ih'v)iy*) - p{xh',Xh' - y*) 

^ {Kh-Uh' - Ah'Ih'v){xh') = (w/i' - Kh'Ih'Vh'){xh') 

and thus 

v[xh') ^ {Kh'Ih'v){xh'). 
The consistency of the discretization, stated in Theorem QiJ yields that 

H{xo,Dv{xo)) < 0, 

which concludes the proof that u is a viscosity subsolution. 

In the same way we prove that w is a viscosity supersolution of H{x, Du{x)) = 0. 
Therefore, u is a viscosity solution, which, by the comparison principle (Theorem [2|, 
is actually unique. Hence, there is exactly one limit point of the sequence Uh, which 
thus has to converge uniformly to the just established viscosity solution u. □ 

Remark. Note that the only use that we have made so far of the existence Theorem Q] 
was to motivate the local variational principle for the finite- element discretization. 
In fact, our proof of the convergence result shows the existence of a viscosity solution 
en route — under the somewhat more restrictive compatibility conditions (H4') and 
(H5'), however. 



7. The Hopf-Lax Update for Generalized Eikonal Equations in 2D 

Let 17 C be a polygonal Lipschitz domain and A/ : ^ M^^^ be a contin- 
uous mapping into the symmetric positive definite 2 x 2-niatrices. We denote the 
corresponding inner product by {p,q)M{x) = {M{x)p,q), its subordinate norm by 

\\p\\m{x) = {p,p)m%)- 

Now, we consider the Dirichlet problem for the generalized eikonal equation, 

\\Du\\m(x) = 1 in fi, u\dn = 9- 

Its Hamiltonian H{x,p) = \\p\\m(x) ^ 1 satisfies the assumptions (H1)-(H3) and 
(H4'). The support function of the zero-level set is simply given by the norm that 
is dual to II • ||M(a;)' namely, 

p{x,q)^ max {p,q) ^ max {p,q) ^ \\q\\M{x)-^- 
The Hopf-Lax update function becomes 

{AhUh){xh) = min [uh{y) + \\xh - y\\M{x^)-n , Xh e Vlh, Uh eVh- 

yedujhixh) ^ ^ 

There is a simple procedure to evaluate (A^u/i)(a;^) at Xh G ^h- To this end let 
fTi, . . . , Um G Tih be the triangles that have Xh as a vertex and Ji the (closed) edge 
of <Ti opposite to Xh- Then, as in (12}, 

{KhUh){xh) = min with ^ miTiiuh^y) + \\xh - y\\M(xh)-A- 
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Figure 2. Geometry ofthe minimization of cos((5)||y—2/h||-|-||a;h—j/|| for 
y & \yh, Zh]- Note that for S > n/2 the segment through y perpendicular 
to Is has negative length cos((5) ■ \\y — yh\\- 



Let us take one of the triangles, ai, (see Figure Q]) and call its vertices Xh, Vh, Zh, 
hence Ji = [yh-, Zh]- In the case of the classic eikonal equation, that is, M{x) = /, 
the update Ui can be determined from an elementary geometric argument. 

Lemma 12. Let a G S/i be the triangle with the vertices Xh, yh, and G Vh- 
Denote the angles at yh, Zh by a, (3, respectively. Defining 

^ _ Uhjzh) - Uhivh) 

\\zh-yh\\ 

and cos((5) — A if \A\ ^ 1, we obtain 

Ui min (uhiy) + \\xh - y\\) Uh{yh) + niin (A 
y&[yh,zh] ^ ' yefeh.zfi] v 

{uh{yh) + \\Xh - yh\\: 
Uh{yh) + cos{S - a) ■ \\xh - yh\\, 
Uh(zh) + \\Xh - Zh\\, 

Proof. For |A| ^ 1 the assertion follows from a direct application of the triangle 
inequality; e.g., for A ^ 1, 

A • ||y - yh\\ + \\xh - y\\ ^ \\y - yh\\ + \\xh - y\\ > \xh - yh^- 

Now, let |A| < 1 so that cos((5) — A defines a (5 G (0,7r). A look at Figure^ shows 
that 

cos(5) • II?/ - + llx,, - y|| (16) 
attains its minimum at the unique intersection of two straight lines: the first line 
running through yh and z/i, the second line running through Xh perpendicular to 1^. 
Here, l(, is the straight line that encloses at yh with \yh^ Zh] the angle Tr/2 — S. We 
observe that the value of the minimum is simply cos((5 — a) ■ \\xh ~ yh\\- A further 
look at Figure \2\ teaches that G [yh , Zh] if and only if 

O^S — a^"f — TT~a~(3, that is, a ^ S ^ tt — p. 

If S < a, or equivalently A > cos(q;), is to the left of yh and the minimum of 
p6l ) in [yh, Zh] is attained at yh- On the other hand, if i5 > tt — /?, or equivalently 
A < cos(7r — is to the right of Zh and the minimum of in [yh,Zh] is 

attained at Zh. □ 



■ \\y ~ VhW + \\xh ~ y\\) 

cos(q;) ^ A, 
A cos(7r - /?). 
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For the general case we simple apply the triangular update formula of Lemma 22] 
to the image of the triangle ai under the linear transformation M{xh)~'^^'^ ■ This 
way we immediately obtain the following update procedure, writing {p,q)x = 
(p,9)M(a;)-i, IblU = lblUf(a;)-i, Ca = cos(a), and Cf3 = cos(/5) for short (note that 
we used the addition formula to spell out cos(q; — S) for implementation purposes): 

^ _ uhjzh) - uh{yh) ^ 
\\zh-yh\\xu 

_ {xh - Uh, Zh - yh)xh _ {xh- Zh^Vh- Zh)xh 

\\xh-yh\\xH-\\zh-yh\\x^' \\xh - zhWxH ■ Wvh - zhWxH 

if Ca A 

Ui = Uh[yh) + \\xh - yhWxh^ 
else if A ^ —cp 

= Uh{zh) + \\xh - ZhWx^] 

else 

= uhiyh) + (c„A + V(l-c2)(l- \\xh - yh\\x^] 

Remark. With different ideas on a discretization, exactly the same update formula 
has been obtained for the (classic) eikonal equation by Kimmel and Sethian [ KS98] 
(see also |Set99i §10.3.1]), who use for acute triangulations the methodology of 
[BS98] to construct upwind schemes on unstructured meshes, and, independently, 
by the geophysicist Fomel [Fom97], who locally uses Format's principle of shortest 
travel times (which is closely related to our local use of the Hopf-Lax formula). 

Sethian [Set99, §10.1] shows further that this update formula generalizes the up- 
wind finite-difference scheme on structured grids given by Rouy and Tourin [RT92] . 



8. Solving the Discrete System 

A Review oj Methods. Theorem [6] shows that the nonlinear discrete system ( ITT] ) can 
be solved by the fixed-point iteration 

ul+^^khK, n = 0,l,2,..., 

for a suitably chosen initial iterate u^. Such a fixed-point iteration uses the updated 
values at a nodal point Xj only after all the updated values have been calculated. 
This corresponds to the classic Jacobi-iteration for linear systems of equations and 
lends itself to direct parallelization. 

If we sequentially traverse the nodal points in a given order and modify the 
iteration to always use the most recently updated value, we obtain a nonlinear 
variant of the Gauss-Seidel iteration. Rouy and Tourin [RT92] used such a non- 
linear Gauss-Seidel iteration to solve a finite difference discretization of the eikonal 
equation ||£'u(a;)|j — n(x) on a structured mesh. 

For both iterative methods the complexity will typically scale as 0(A^^+-'^/''), 
where N denotes the number of nodal points and d the space dimension. This 
is because the information about the solution, inherent initially to the boundary 
only, travels by next neighbor interaction at each run trough all nodal points. To 
spread that information to the whole computational domain about 0{N^^'^) runs 
are necessary. However, even though this heuristic well explains the experimental 
observations, to our knowledge there is no rigorous proof of that in the literature. 
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In 1995 Sethian [Set96] and Tsitsiklis [ Tsi9 5j have shown independently that 
for eikonal equations on structured meshes the nonhnear equation can be solved 
exactly in a single pass, that is, by traversing the grid once using local operations 
only. This fast marching method was later generalized to triangular meshes by 
Kimmel and Sethian [KS98]. It relies on the causality property, namely that on 
an acute triangulation the value Uh{xh) depends only on the values in neighboring 
nodal points yu that are lower, Uh{yh) ^ Uh{xh)- So the discrete solution can be 
computed starting from the neighborhood of the boundary moving further inwards 
the computational domain along increasing values of Uh- However, on non-acute 
triangulations additional effort is necessary to deal with the loss of this causality 
property (see [KS98] for details). The complexity of this method is 0{N \og{N)) 
where the logarithm comes from administering a priority queue of candidates for 
the next lowest value of Uh, such as a heap data structure. 

For the particular Hamilton-Jacobi-Bellman equation (8]) a single pass algo- 
rithm generalizing the fast marching method, called the ordered upwind method 
(OUM), was introduced by Sethian and Vladimirsky [SVOOJ and is discussed in 
detail in [SV03 ] . This method is not a fast solver for a given discretization, but 
the discretization is specifically designed for the needs of the fast solver. Like the 
Hopf-Lax update the update formulas of the OUM are based on local variational 
principles (Bellman's principle). However though, the OUM does not solve (|TT|). 
since the update in Xh is not necessarily computed from the neighborhood ujh{xh) 
but from larger neighborhoods within the radius h ■ v. Here v = F*/F^, denotes 
the anisotropy coefficient of the Hamilton-Jacobi-Bellman equation (see (7])). This 
quantity v affects not only the complexity of the OUM, which is 0{v'^^^ N \og{N)) , 
but also its accuracy. 

Adaptive Nonlinear Gauss-Seidel Iteration. In this paper we propose an adaptive 
version of the nonlinear Gauss-Seidel iteration, which is a modelled after a similar 
relaxation method [ PR93] for the multilevel solution of elliptic boundary value 
problems. It turns out to be substantially faster than the standard Gauss-Seidel 
iteration, easy to implement and universal. 

The adaptive Gauss-Seidel iteration differs from the standard one in two respects. 
First, like in the fast marching method, only those nodal points are updated that 
"have the information", that is, are neighbors of recently updated points. Second, 
the order of updates is not fixed but varies as the iteration proceeds. Thus, a queue 
denoted by Q is administered to provide the ordering of updates. However, other 
than in the fast-marching method where using the causality property requires to 
keep control of the point with minimal function value, the queue is now simply 
given the structure of a FIFO (first in first out) stack: the nodal point staying 
longest in the queue is updated next. 

The algorithm is passed a user-defined tolerance tol and it ends up with an 
approximate finite-element solution Uh G Vh such that 

\\uh ~ KhUhW^ ^ tol. 

It is organized as follows: 

(1) (Initialization) Let M/i|an^ = fflaOh, i*ft,|nh = '^■^ Let Q be the list of all 
points Xh £ that are adjacent to some boundary point (in an arbitrary 
but fixed order). 

(2) (Iteration) Remove the first point x^ from Q and compute the update value 

Unew = {AhUh){Xh)- 



In fact any value larger than the bound ( |14)) will do. However, taking oo makes the argument 
more elegant and can correctly be implemented in IEEE arithmetic. 
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(3) If |uncw — Uh{xh)\ > tol then update Uh{xh) = Wnew and append all not yet 
enqueued neighbors yh of Xh to the queue Q. 

(4) If Q 7^ 0, goto (2). 

To prove the convergence of this method, we denote the initial finite-element 

function of step 1 by After the nth update has been performed in step 3 the 
actual finite-clement function will be denoted by u^. 

Theorem 13. The algorithm generates a sequence u^, u\, . . ., that is monotonously 
decreasing. It terminates after finitely many steps with an approximate finite ele- 
ment solution Uh € Vh, such that \\uh — ^hUh\\ ^ tol. 

Proof. The initialization Uhln^ = oo ensures that every point Xh is updated at least 
once, as the residual is oo when the first update value in Xh is computed. After 
the first update, Uh{xh) is assigned a finite value, since Xh has a neighbor in d^h. 
or a neighbor, for which a finitely valued update has already been computed. By 
induction on n we get that at each later update of a nodal point Xh, all neighbors 
of Xh that have been changed over the last update can only have been assigned a 
lower value of Uh- ;,From the monotonicity of Kh we thus get the first assertion. 

Since an update in step (3) only affects the residual in the neighboring points, 
which are are immediately enqueued, it holds that 

{xhe0.h : < oo and |AftU^ - u^l > tol} C Q 

for every n ^ 0. So if the algorithm terminates with Q = 0, the tolerance has been 
reached. 

Otherwise, if the iteration does not terminate, then there is at least one nodal 
point xjjj that appears infinitely often as the first element of the queue Q and gets 
updated at steps Uj — > oo, j — > oo. Hence, there must be \u^' {^h)~''Hi'~^{^h)\ > 
in contradiction to the convergence of m)^^ (x^) as j ^ oo which is implied by the 
monotonicity and the trivial lower bound ^ vamx^ao, g{x). □ 

Though the run-time complexity of the adaptive Gauss Seidel iteration behaves 
probably at worst as 0(7V^+^/'') like in the standard Gauss-Seidel iteration, a lot 
of unnecessary updates are saved as we will see in the numerical experiments of the 
next section. 
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Figure 4. Left: Value function of some min-time optimal control prob- 
lem. Right: Complexity /accuracy of the methods in comparison. 



9. Numerical Experiments 

For the two following examples the solutions were computed on unstructured 
iTieshes with 23^, 45^, 91^, 181^, and 725^ nodal points. A solution on a mesh with 
1451^ points served to estimate the discretization error. The iterative methods were 
used with an absolute tolerance tol = 10^*. 

The first example concerns the distance map on the torus given by the immersion 

f{xi,X2) = (cos(27ra;i)(5 + 4cos(27rx2)), sin(27rxi)(5 + 4 cos(27ra;2)), sin(27rx2)) . 

Using the Gram matrix G{x) = Df{x)^Df{x) and p'^{x,q) = {q,G{x)q) the dis- 
tance between the points f{x) and /(j/) on the manifold is given by the function 
6{x, y) as defined in |2). With the results of §[7] and }Lio82[ Thm. 5.3(iv)] we obtain 
that u{x) = 5(x,0) is the viscosity solution of the Dirichlet problem 

=lonr!\{0}, u(0) = 0, 

where D, = [—0.5,0.5]^. The solution is shown to the left of Figure S] as a contour 
plot. 

To the right of Figure [3] we compare the accuracy and complexity of the adaptive 
Gauss-Seidel method with both the standard Gauss-Seidel iteration and the GUM.'* 
Here, by complexity we mean the total number of updates calculated on a triangle 
by a formula such as the one at the end of JZI We observe that the adaptive Gauss- 
Seidel iteration is more than a factor of 10 faster than the standard Gauss-Seidel 
iteration but displays the same asymptotic rate of complexity. The GUM show, as 
theoretically expected, a better rate of complexity that, in this example, becomes 
significant even at larger tolerances. 

The second example is taken from [SVOl] and shows the effect of a moderately 
large anisotropy coefRcient v = 19, which affects the complexity of the GUM. We 
consider a simple min-time optimal control problem governed by the dynamical 
system 

y'{t)^a{t) + h{y{t)), y{Q) = x. 
The controls a(-) are taken from A~ {a: [0,oo) measurable} and 

y 

b{y) — — 0.9 sin(47r2/i) sin(47r?/2) ■ -jr-jr- 



We have coded the OUM from [SV03] with a little completion that turned out to be necessary: 
Considered points have also to be updated, if they depend on an edge that drops out of the accepted 
front. If some point Xh gets accepted this may happen to any edge opposite to xj^ in u)fi{xf^). 
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For X G fl — [—0.5, 0.5]^ and a control a we denote by Tx{a) the minimal time that 
the trajectory y{-) takes to reach the origin. 

Following [BCD97i p. 241, Thm. 2.6] the value function u{x) = MaeATxia) is 
the viscosity solution of the Hamilton-Jacobi-Bellman equation 

H{x, Du) ^ max (a + b{x),-Du) -1 = 0, it(0) = 0. 

||a|| — 1 

One figures out that H{x,p) = ||p|| - {b{x),p) - 1; (HI), (H2), (H4'), and, as 
||5|| ^ 0.9, the coercivity condition (H3) are fuUfilled. A short calculation shows 
that 



[l-\\bix)f + {b{x),q/\\q\\fy^'-{b{x),q/\m 

compare also (S) and [SVOll Eq. (19)]. The solution calculated on a 253 x 253 
mesh can be found to the left of Figure 31 To the right of this figure the accuracy 
of the approximate finite-element solution is shown versus the complexity of the 
iteration. A comparison of the (adaptive) Gauss-Seidel iteration with the OUM is 
shown to the right of Figure SI Again we observe that the adaptive variant of the 
Gauss-Seidel iteration is by about a factor of ten more efficient than the standard 
one. This time, however, the quite sophisticated order upwind method behaves less 
favorable: because of the large anisotropy coefficient the break-even point at which 
the OUM becomes more efficient than the simple adaptive Gauss-Seidel iteration is 
at a mesh-size of more than 725^ = 525,625 nodal points. We expect this effect to 
become even more pronounced in 3D and higher, because of the increasingly better 
complexity rate of the Gauss-Seidel iteration. 
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